﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.ComponentModel.Composition;

namespace Nature.Numerics
{
    /// <summary>
    /// The AHMms2004 implementation of the minimum mean square error algorithm 
    /// </summary>
    [Export(typeof(IFuncBasisResolver))]
    public sealed class AHMms2004FuncBasisResolver : IFuncBasisResolver
    {
        /// <summary>
        /// Initializes a new instance of the AHMms2004FuncBasisResolver class.
        /// </summary>
        public AHMms2004FuncBasisResolver()
        {
            this.Solver = new AHRpGauss2004();
        }

        /// <summary>
        /// Initializes a new instance of the AHMms2004FuncBasisResolver class.
        /// </summary>
        /// <param name="solver">Linear algebraic equations solver used in solving the function extremum problem</param>
        [ImportingConstructor]
        public AHMms2004FuncBasisResolver(
            [Import(typeof(ILinearAlgebraicEquationsSolver), AllowDefault=true)]
            ILinearAlgebraicEquationsSolver solver = null)
        {
            this.Solver = solver ?? new AHRpGauss2004();
        }
        
        /// <summary>
        /// Gets the instance of the linear algebraic equations solver
        /// </summary>
        public ILinearAlgebraicEquationsSolver Solver { get; private set; }

        /// <summary>
        /// Finds the functional basis decomposition coefficients
        /// </summary>
        /// <param name="function">The function to be represented in the form of the functional basis</param>
        /// <param name="domainPoints">Collection of the independent variable values (x-domain)</param>
        /// <param name="basisCollection">Collection of functions representing the functional basis</param>
        /// <returns>Linear decomposition of the given function</returns>
        public FuncBasisDecomposition Decompose(
            Func<double, double> function, 
            IEnumerable<double> domainPoints,
            IEnumerable<Func<double, double>> basisCollection)
        {
            double[] domain = domainPoints.ToArray();
            Func<double, double>[] basis = basisCollection.ToArray();
            int nva = basis.Length;
            int npo = domain.Length;
            double[][] a = NumericalTypeConverter.AllocateMatrix(basis.Length);
            double[] b = new double[nva];
            for (int i = 0; i < nva; ++i)
            {
                double sum_b = 0.0;
                for (int ipo = 0; ipo < npo; ++ipo)
                {
                    double xk = domain[ipo];
                    double fk = function(xk);
                    double fi = basis[i](xk);
                    sum_b += fk * fi;
                }
                b[i] = sum_b;

                double[] row = a[i];
                for (int j = 0; j < nva; ++j)
                {
                    double sum_a = 0.0;
                    for (int ipo = 0; ipo < npo; ++ipo)
                    {
                        double xk = domain[ipo];
                        double fk = basis[j](xk);
                        double fi = basis[i](xk);
                        sum_a += fk * fi;
                    }
                    row[j] = sum_a;
                }
            }
            Solver.Solve(a, b);
            return new FuncBasisDecomposition(basis, b);
        }
    }
}
